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Abstract — We analyze the behavior of the outer envelope in a massive star during and after the collapse 
of its iron core into a protoneutron star (PNS) in terms of the equations of one-dimensional spherically 
symmetric ideal hydrodynamics. The profiles obtained in the studies of the evolution of massive stars up to 
the final stages of their existence, immediately before a supernova explosion (Boyes et al. 1999), are used as 
the initial data for the distribution of thermodynamic quantities in the envelope. We use a complex equation 
of state for matter with allowances made for arbitrary electron degeneracy and relativity, the appearance 
of electron— positron pairs, the presence of radiation, and the possibility of iron nuclei dissociating into 
free nucleons and helium nuclei. We performed calculations with the help of a numerical scheme based 
on Godunov's method. These calculations allowed us to ascertain whether the emersion of the outer 
envelope in a massive star is possible through the following two mechanisms: first, the decrease in the 
gravitational mass of the central PNS through neutrino-signal emission and, second, the effect of hot 
nucleon bubbles, which are most likely formed in the PNS corona, on the envelope emersion. We show 
that the second mechanism is highly efficient in the range of acceptable masses of the nucleon bubbles 
(<O.O1M0) simulated in our hydrodynamic calculations in a rough, spherically symmetric approximation. 
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INTRODUCTION 

The root cause of supernova explosions in massive 
stars (Mjvis > IOM0) after the completion of their 
thermonuclear evolution is commonly assumed to be 
the gravitational collapse of their central iron cores 
with masses Mpe in a narrow range (I.2M0 < Mpe < 
2Mq); the lower limit is approximately equal to the 
Chandrasekhar mass of an iron white dwarf. To be 
more precise, iron-core collapse is followed by a num- 
ber of hydrodynamic processes in the compact central 
part of the cavity left from the iron core. 

Hydrodynamic spherically symmetric models for 
the gravitational collapse of iron stellar cores [see, 
e.g., the review article by Imshennik and Nadyo- 
zhin (1982)] showed a clear separation of the collapse 
into two stages. At the first stage, the inner iron 
core with mass Mjpe < lAf© (Imshennik 1992) or, 
to be more precise, with mass (0.6— O.8)M0 (Nady- 
ozhin 1998) collapses homologically. The remain- 
ing outer iron core with mass MePe = Afpe — Mjpe — 
(0.4— 1.4)Mq in the above range of iron-core masses 
lags well behind the inner core in contraction param- 
eters (density, etc.). In a rough approximation, it can 
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even be assumed to maintain the hydrostatic equilib- 
rium of the initial state. Since these issues were dis- 
cussed previously (Imshennik and Zabrodina 1999; 
Imshennik and Popov 2001), the hydrodynamic be- 
havior of the outer iron core and outer layers of a 
massive star may be called the problem of their post- 
shock accretion (Brown et al. 1992) onto the embryo 
of a protoneutron star (PNS) with a mass of about 
IMq formed at the end of the first collapse stage. The 
second collapse stage of the matter surrounding this 
PNS embryo can be identified with this post-shock 
accretion. Literally, this implies that at the initial time, 
there is a collapsed central object with mass IM© 
and its surrounding layers in an initially hydrostatic 
equilibrium. These initial conditions for the model of 
post-shock accretion considered below seem to be 
in qualitative agreement with hydrodynamic models 
for the collapse of iron stellar cores. They allow the 
hydrodynamic processes in the remaining onionlike 
part of the star surrounding the iron core in a massive 
star of arbitrary total mass (M^s > IOM0) to be 
included in the analysis. Here, we show that these 
processes do not result in an explosion of the outer 
stellar layers surrounding the iron core as during a 
supernova explosion (~10^^ erg). However, in the 
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presence of low-mass (<O.O1M0) nucleon bubbles, 
they cause the emersion of the outer layers with a 
relatively low energy and with a total mass of several 
Mq up to the inner boundary of the helium shell. 

FORMULATION OF THE PROBLEM 

The System of Equations 

In most cases, the hydrodynamic behavior of 
the envelope in a massive star can be described 
by the system of equations of ideal hydrodynam- 
ics. In our spherically symmetric case 

d d \ 

^,50 = 0,5<^ = 1 , this system m spherical 

coordinates (r, 9, (p) is known to be 



dp 19,9^^, 
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at ■V^d-yp''r^^^ = P^- 

^ + ^2§-yyr{pE + P))=pVr9r. (3) 

Here, E = e + — is the sum of the specific inner 

(s) and kinetic energies. The radial acceleration 
of gravity gr is the sum of two components: gr = 
9pns + 5env The first component is attributable to 
the gravitational field produced by the PNS located 
exactly at the coordinate origin: 
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The second component is attributable to the intrinsic 
gravitational field of the stellar envelope. The acceler- 
ation of this field is defined by the standard equation 
genv = — V$ and the potential satisfies the Poisson 
equation 

A$ = AirGp. (5) 

In our calculations, we actually used a numerical 
scheme based on the equations of ideal hydrodynam- 
ics written in axisymmetric form ( g^, = ] in 

\dip ) 

spherical coordinates (r, d, tp): 
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As can be easily verified, this system transforms to the 
system of equations (1)— (3) if we set Vq, V^p, and ge 
equal to zero. 

Equations (6)— ( 10) are written in divergence form; 
since they do not contain viscosity, including artifi- 
cial viscosity, shock and contact discontinuities can 
emerge in their solutions. The absence of heat con- 
ductivity in these equations makes discontinuities in 
specific internal energy and density admissible. This 
can be taken into account when choosing a finite- 
difference scheme. 

Initial Data 

As we noted in the Introduction, the formation 
time of the PNS embryo with the mass {^IMq) 
characteristic of our problem, when the neutron sig- 
nal also increases most steeply (even before its max- 
imum), may by arbitrarily taken as the initial time 
{t = 0). In that case, the outer layers may still be 
considered in a hydrostatic equilibrium. We used the 
distributions of thermodynamic quantities obtained in 
the studies of the evolution of massive stars (Boyes 
et al. 1999) as the initial data. From this paper, we 
took the density and temperature profiles, which, in 
turn, serve to determine the initial pressure and inter- 
nal energy by using a specified equation of state. Of 
course, the equation of state used in our calculations 
slightly differs from the equation of state used in ob- 
taining these profiles. The main difference is that we 
took a large number of isotopes into account in our 
numerical calculations of stellar evolution. However, 
a detailed comparison of the equations of state shows 
that, quantitatively, these differences are moderately 
large and passing to a simpler equation of state in 
our calculations has no significant effect on the radial 
distributions of thermodynamic quantities. The pres- 
sure calculated from fixed density and temperature 
using our equation of state differs by no more than 2% 
along the entire profile from the pressure taken from 
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Boyes et al. (1999); the pressure from Boyes et al. 
(1999) is lower approximately by 2% near the iron 
core and higher by the same ~2% in the profile tail. 
Nevertheless, to exactly specify hydrostatically equi- 
librium distributions of thermodynamic quantities in 
our calculations as the initial data, we recalculated the 
initial profiles using a new equation of state. As was 
noted above, by the formation of a PNS embryo with a 
typical mass of ^IMq, the outer layers surrounding 
the iron core are virtually in hydrostatic equilibrium, 
as suggested by hydrodynamic collapse calculations 
(Nadyozhin 1977). To be more precise, the radial 
temperature distribution from Boyes et al. ( 1999) was 
used in these calculations and the new distributions 
of pressure and the remaining thermodynamic quan- 
tities (p, e) were reconstructed by solving the system 
of hydrostatic equilibrium equations 

(11) 



dr ' 

dm 2 
— = 47rr p- 



(12) 



the pressure P and mass m at the inner boundary 
were also taken from Boyes et al. ( 1 999). The dimen- 
sions of the computed region were chosen in such a 
way that the inner boundary was equal to the radius 
that bounded the region of mass IMq in the initial 
profile and the outer (in radius) boundary coincided 
with the inner boundary of the helium shell. The cal- 
culations were performed for the evolutionary models 
of stars with masses of 1 1, 20, and 25M0 and solar 
metallicity Z. 

The initial values of the velocity component Vr 
were taken to be zero in all cases, except for an 
additional calculation in which we imparted a velocity 
approximately equal to the free-fall velocity in magni- 
tude and directed away from the center. 

In all of our calculations, we took into account 
the presence of a neutrino signal. The gravitational 
collapse of massive stellar cores gives rise to a 
powerful neutrino signal, as predicted by Nady- 
ozhin (1978) and experimentally confirmed by the 
observations of the SN 1987A explosion in the Large 
Magellanic Cloud. The neutrino signal carries away 
almost all of the gravitational binding energy of the 
forming neutron star, which is released during the 
collapse. The total gravitational mass defect of the 
neutron star was calculated by using the formula 

AMg = 0.084-^ (Lattimer and Yahil 1989), which 

is AMg = 0.054, 0.165, and O.SSOM© for iron cores 
with masses Mpe = 0.8, 1.4, and 2.OM0, respectively. 
It should be immediately noted that these values of 
AMg refer to cold neutron stars; i.e., they definitely 
overestimate the effect under consideration, because 



a hot PNS still has high internal energy and tem- 
perature. The shape of the neutrino-signal cur\'e in 
time was taken from Nadyozhin (1978). During the 
calculation, it was fitted by a set of exponentials that 
were continuously joined at the boundary points. 

Boundary Conditions 

The region of the solution of the problem or the 
computed region has the shape of a spherical envelope 
(see above), rmin <r< rmax; the choice of rmin and 
''max is determined by the physical considerations 
outlined in the Subsection "Initial Data." Sufficient 
boundary conditions must be specified precisely on 
these constant values of the Eulerian radius. 

Thus, in our problem, the computed region has 
the Eulerian outer and inner boundaries, rmin and 
''max- The inner boundary is a transparent wall; i.e., 
the derivatives of all physical quantities (vr,p, and 
e) are assumed to be zero. Matter can then freely 
flow through the inner boundary. In that case, the 
accretion rate or, more precisely, the total mass of the 
matter passing through it per unit time is calculated 
at this boundary. At the outer boundary, the condition 
simulates a vacuum outside the computed region; i.e., 
the thermodynamic quantities {p,e, and P) are set 
equal to nearly zero. The boundary condition for the 
gravitational potential is $ ^ —GM jr for r ^ 00. 
Recall that we use the numerical scheme of a two- 
dimensional problem (see the Section "The System 
of Equations"). 

THE EQUATION OF STATE 

We use the equation of state for matter treated 
as a mixture of a perfect Boltzmann gas of nuclei, 
a perfect Fermi— Dirac electron— positron gas, and 
blackbody radiation. In a wide temperature range, the 
matter is assumed to be a mixture of a perfect iron 
gas of nuclei, a perfect electron— positron gas, and 
blackbody radiation. However, in order that the equa- 
tion of state be applicable under nuclear-statistical- 
equilibrium (NSE) conditions, i.e., at high tempera- 
tures T > Tc, where Tc = (3-5) x 10^ K(lmshennik 
and Nadyozhin 1965, 1982), virtually independent of 
the matter density p (Imshennik et al. 1981 ), we as- 
sume that at T > Tc, the matter is a mixture of perfect 
gases of four nuclides, Jn, \p, aHe, leFe, with a perfect 
gas of e~, e+ leptons and blackbody radiation. In- 
cluding many other nuclides should not significantly 
affect the thermodynamic functions of the matter. 

The contribution of the electron— positron com- 
ponent is given by general integral equations with 
arbitrary degeneracy and relativity. The difference be- 
tween the electron and positron number densities 
can be determined from the necessary condition of 
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electrical neutrality. In turn, the latter in the form 
of a simple integral equation serves to calculate the 
electron chemical potential (//): 



oo 



(13) 



_P_ 

m,. 



where Xpe, -^He> and Xp are the mass fraction of 
iron, helium, and free protons; and F± (^) under the 
integral are the Fermi— Dirac functions. For arbitrary 
degeneracy and relativity of 6=*= leptons, these func- 
tions can be expressed explicitly: 

i^±(e) = e'(i + expz± 



(14) 



because ii = ii- 



-fi+. Here, we introduced the di- 

mensionless parameters a = and w = ^ „ , 

k^T meC^ 

where T is the temperature, is the Boltzmann 

constant, is the electron mass, is the chemical 

potential (including the fermion rest energy), and the 

quantity A as the natural unit of length (of the order of 

3 
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h 



the Compton length) = — , 

The equation of state proper, i.e., the expressions 
for the pressure and specific internal energy of the 
matter, is then (Imshennik and Nadyozhin 1965, 
1982) 



P = P- + P+ + \aT^ 

H —pT (Xn + Xp+ -Xue + T^Xpe 

ruu \ 4 56 



(15) 



aT'^ 3 % _ 

= £_+£+ + + -— T 

p 2mu 

Xn + Xp + ^-^He + 5g"^Fe 



(16) 



+ 



Qpe + 26 AQ^ 
56m„ 



(1 - Xpe) 



QHe + 2AQ„ AQn ^ 
-AHe -^n 



4m,, 



mp) = 1.294 MeV is the en- 



Here, AQn = {rrin 

ergy threshold of P decay and a = 7.5644 
10"-^^ erg cm~^ is the neutron radiation density 
constant. 

The contributions of the lepton pressure and spe- 
cific internal energy to the complete equation of state 



are given by standard general relations (Landau and 
Lifshitz 1976) [with definitions (14)]: 



e 



rrieC 

^ 



[F-iO + F+iOH, 

(17) 



meC 





(18) 



UIqC 



26 1 
56 m. 



Note that a constant electron rest energy at zero tem- 
perature (in the absence of positrons) is subtracted 
from the sum of e_|_ and £_. Clearly, the specific in- 
ternal energies can generally be determined to within 

an arbitrary constant. 

The mass fractions of iron, helium, free neutrons, 
and free protons under NSE conditions at a given 
density and temperature are defined by the system 
of equations (Imshennik and Nadyozhin 1965, 1982; 
Imshennik and Zabrodina 1999) 



Xi = LOiA^/^ 



P 

rUu 



2'Kmuk^T 

Ai-l 



(19) 



. Xn + (l/2)XHe + (30/56)XFe 

t'O = -TP , /-, /ON ^ , /r,nlrnw = 



30 



Xp + (l/2)XHe + (26/56)XFe 26 

Xn + X„ + + ^Fe = 1, 



(21) 



where i = He, Fe; m„ is the atomic mass unit; Qi is 
the binding energy in the ground state of the nuclides 
(QHe = 28.296 MeV, Qpe = 492.26 MeV); Ai and Zi 
are the mass and charge numbers of these nuclides; 
and uji are the partition functions, which may be taken 
to be unity without any particular error. Here, we as- 
sumed the ratio of the total numbers of neutrons and 
protons to be constant; 9 = 30/26 = ^pe is the value 
typical for the IgFe nuclide (Imshennik and Zabrod- 
ina 1999). Actually, this implies that we completely 
ignore the (3 processes and disregard the variety of 
iron-peak elements among the explosion products of 
a low-mass neutron star and in the surrounding iron 
gas. Clearly, only iron nuclei are available at T < Tq; 
i.e., Xpe = 1 and X^^ = Xn = Xp = Q. 

The term in the square brackets from (16) is the 
rest energy of the matter as a function of the mass 
fractions, which, obviously, becomes zero at Xpe = 1. 
This means that the energy of pure iron is taken as 
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a zero level, being, of course, the lowest value of the 
above term. 

Since the pressure and the specific internal en- 
ergy must regularly be used in solving the system of 
hydrodynamic equations, it would be appropriate to 
use a tabulated equation of state. This would allow 
us to avoid complex calculations at each step of the 
solution of system (6)— (10) if the thermodynamic 
functions were calculated once and with a high ac- 
curacy. Our numerical method for solving the system 
of hydrodynamic equations (6)— (10) uses the density 
and specific internal energy as the main variables. 
Therefore, the equation of state is P = f{p,e). In 
the numerical method (see below), this equation is 
simulated by the so-called binomial equation of state 



P=[{^-l)e + cl]p- pocg. 



(22) 



This simulation is local in nature. The constants 7, 
Cq, and po are determined by the pressure and its 
derivatives at a given point; 

-^(^) 



(23) 



POCo 



, df 

pya-p 



f, 



1 

7-1 = - 

p 



dl 
de 



By directly substituting (23) into (22), it can be 
verified that the simulation of equation of state (22) 
with specified [according to (23)] constants satisfies 
the requirement that the pressure and the speed of 
sound derived from (22) be equal to those derived 
from the real equation of state P = f {p,e). Using 
the binomial approximation for the equation of state 
considerably simplifies the solution of the problem of 
discontinuity breakup, which underlies the numeri- 
cal method. Thus, for the system of hydrodynamic 
equations to be numerically solved, we must know 
the pressure and its first derivatives with respect to 
density and specific internal energy as a function of 
these variables. 

Equations (15)— (16) do not allow us to obtain P 
as a function of p and e in explicit form. Therefore, it is 
necessary to first determine all of the thermodynamic 
quantities concerned as functions of density and tem- 
perature. To this end, we used an approach related 
to the Gaussian quadrature method (Krylov 1967). 
This approach has a sufficient accuracy over the en- 
tire range of densities and temperatures concerned 
(Blinnikov et al. 1996). To solve the system of hy- 
drodynamic equations (6)— (10), we must also know 
the first derivatives of the pressure with respect to 
density p and specific internal energy e. Using stan- 
dard relations for the derivatives of thermodynamic 
quantities (Landau and Lifshitz 1976), we find the 



thermodynamic derivatives of the pressure needed for 
the numerical scheme to be implemented 



dP\ 





^ > (24) 



de\ \dp J rp 



dP\ 




where P{p,T)=P{p,e{p,T)). 

Finally, for the equation of state to be derived in 
final form, we must reinterpolate the computed tables 
of thermodynamic quantities and their derivatives, 
which depend on density and temperature, to the 
domain of the p and e variations. 

Clearly, using our equation of state, we can sim- 
ulate the formation of a nucleon bubble within the 
computed region by artificially raising the tempera- 
ture in the initial conditions compared to the temper- 
ature obtained by Boyes et al. (1999). In the above re- 
calculations of the hydrostatic equilibrium conditions 
with the inevitable decrease in density on the sim- 
ulation segments, we need not be concerned about 
the conservation of a purely nucleon composition, 
because the decrease in density only facilitates it. 

THE NUMERICAL METHOD OF SOLUTION 
The numerical solution of our problem is based on 
the popular and universal PPM method. This method 
uses the Eulerian finite-difference scheme (Colella 
and Woodward 1984) and is a modification of Go- 
dunov's method. 

The numerical calculation was performed in spher- 
ical coordinates. To construct the finite-difference 
scheme, we rewrote the initial system of equa- 



in 



tions (6)-(10) 
{p,pVr,pVe,pV^,pE}'^: 

nF,(U)) 



dV d{Ar 
'dt^ dDr 
dG 
dr 

( 



variables 



9(^eFe(U)) 
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pVrVe 
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\pVeE + VgPj 



(25) 
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\V 
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J = 



^ 

P9r'rp{Vi + V^)lT 

pge + p{V^coW - VrVe)/r 
-pV^iVr + VecoW)/r 
p{grVr + geVe) j 



Ae = 



Dg = — COS 9. 



System (25) is hyperbolic; to solve it, we broke down 
the computed region into cells using a grid composed 
of surfaces of constant radius (r) and constant po- 
lar angle (6). Next, we formally averaged the equa- 
tions of system (25) over the cell volume AVijk and 
over the time step At = t""*"^ — t". The problem of 
determining the distribution of the physical quanti- 
ties {p, pVr, pVg, pV^, pE}ij at a later time t^^^ from 
their known distribution at a given time is then 
reduced to determining the time-averaged fluxes at 
the cell boundaries and the time- and cell-volume- 
averaged free terms if, of course, the latter exist. 
To calculate the average fluxes, the interaction be- 
tween two adjacent cells is considered as the prob- 
lem of discontinuity breakup, whose solution is the 
main procedure of our numerical method. In gen- 
eral, an arbitrary discontinuity self- similarly breaks 
up into a configuration composed of four flow regions, 
which differ in properties and which are separated 
from one another by shock (or rarefaction) waves 
and a contact discontinuity (see, e.g.. Landau and 
Lifshitz 1986; Kibel et al. 1963). In constructing the 



finite-difference scheme, we simulated the equation of 
state by the binomial formula (22). The use of this 
formula made it possible to properly describe strong 
rarefaction waves, which is of particular importance in 
astrophysical problems characterized by density vari- 
ations over a wide range and by large density jumps. 
In this case, we can avoid numerically solving a sys- 
tem of ordinary differential equations and, without re- 
sorting to the formal replacement of rarefaction waves 
by shock waves (Colella and Glas 1985), can easily 
write out analytic formulas for shock and rarefaction 
waves. These formulas are required to construct an 
iterative process of calculating the pressure and ve- 
locity in the central region formed after a discontinuity 
breakup. 

Godunov's method can be generalized to a mul- 
tidimensional case (Godunov et al. 1976). However, 
the difficulty in determining the regions of influence 
for the problem of discontinuity breakup arises for our 
algorithm based on the Eulerian scheme. A possible 
solution of this problem is that several regions of 
influence, each corresponding to a certain character- 
istic (Riemann invariant), are found for the boundary 
between adjacent cells. The interpolation distribu- 
tions of the physical quantities are averaged over all 
the derived regions of influence and the initial data for 
the problem of discontinuity breakup are obtained as 
the sums of these means with weighting coefficients 
proportional to the absolute values of the increments 
of the Riemann invariants along the corresponding 
directions. As a result, the characteristic properties 
of the system are properly taken into account (Dai 
and Woodward 1997). We use a parabolic interpola- 
tion for the physical quantities within a cell, which 
is monotonic and continuous at the cell boundaries 
in the region of smooth flow and conserves the to- 
tal volume integrals for the quantities being interpo- 
lated (Colella and Woodward 1984). The constructed 
scheme is explicit; there is a restriction on the time 
step, the Courant condition, to maintain the stabil- 
ity of the calculation. In addition, the method does 
not require a special choice of artificial viscosity for 
each specific calculation. A finite-difference discrete 
approximation of the functions that describe the fields 
of the physical quantities automatically smears and 
smoothes out the discontinuities. 

Below, we also give the system of units used in our 
calculations. The scales of the physical quantities in 
this system are 



[r] = i?o, [Vr] = m = [V^] = {GMo/Ih 



,1/2 



(26) 



[p] = Mo/(47ri?3), [t] = rT/{GMo 



vl/2 



[P] = GM^/{AnRf 



[E] = GMo/Ro, [T] = (GMoV(47ri?^a,))V4, 
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where Rq and Mq are some length and mass scales. 
In units (26), the hydrodynamic equations with New- 
tonian gravitation contain no dimensionless param- 
eters. It would be natural to use the following char- 
acteristic values from our problem as Rq and Mq: 
Ro = 10^ cm and Mq = 10^^ g. The numerical values 
for the scales of the physical quantities from (26) are 
then 

[r] = 10^ cm, (27) 
[Vr] = [Ve] = [V^] = 2.583 x 10^ cm s"\ 
[p] = 7.958 X 10^ gcm"^ [t] = 3.871 x 10"^ s, 
[P] = 5.310 x 10^^ ergcm-^, 
[E] = 6.674 X 10^^ erg, [T] = 2.894 x 10^ K- 

An efficient algorithm is used to solve the Pois- 
son equation (5) and to determine the gravitational 
acceleration g . This algorithm is convenient to use in 
the finite-difference scheme for integrating the hydro- 
dynamic equations on a stationary grid in spherical 
coordinates (Aksenov 1999). The method is based on 
the expansion of the potential in integral representa- 
tion, 

$ = -G / PfAdr', (28) 
J |r'-r| 

by defining the Legendre polynomials in terms of 
the generating function. Using the addition theorem 
for associated Legendre polynomials, expression (28) 
can be reduced to a form convenient for averaging 
over the cell volume in spherical coordinates. This 



procedure can be used in problems with any number 
of measurements and it is efficient: it requires ~imax 
operations per cell, where Zmax is the number of asso- 
ciated Legendre polynomials used in the calculation. 
In addition, the boundary condition $ — > —GM/r for 
r — >■ oo is automatically satisfied for the potential. 

NUMERICAL RESULTS 

Let us turn to the results obtained when numer- 
ically integrating the system of hydrodynamic equa- 
tions (6)— ( 1 0) with the initial data and boundary con- 
ditions given above. We performed two series of cal- 
culations: in the first series, we only took into account 
the effect of the decrease in the central gravitational 
mass due to the emission of a neutrino signal; in the 
second series, apart from allowance for this effect, we 
also analyzed the influence of hot nucleon bubbles 
that are presumably formed in the PNS corona. 

In the first series of our calculations, we considered 
the following models of presupernova stars: [1] the 
total stellar mass Mt^t = ll-M©, the iron-core mass 
(the mass of the central stellar region in which the 
mass fraction of iron-peak nuclei is more than 80%) 
Afcore = 1.266M0, the gravitational iron-core mass 
defect AMg = 2.676 x lO^^ ^he mass of the matter 
inside the computed region Mjn = O.627M0, the 
inner radius of the computed region rmin = 6.250 x 
10^ cm, its outer radius rmax = 1-937 x 10^ cm, 
and the mass of the central region (up to the ra- 
dius r^ij Meenter = O.9987M0; [2] Mtot = 2OM0, 
Mcore = 1.485Mo, AMg = 3.552 x lO^^ g, = 
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2.825M0, rmin = 6.952 x 10^ cm, r^ax = 6.308 x 
10*^ cm, and Mcenter = O.9992M0; [3] Mtot = 25M0, 

Mcore = I.62OM0, AMg = 4.386 X 10^2 g, Min = 

5.458M0, rmin = 8.764 x 10^ cm, r^ax = 6.607 x 
10^ cm, and Mcenter = I.OOOIM0. In all of these 
calculations, we took the constant, relatively low 
specific internal energy £5 = 10^^ erg g"-*^ and density 
Pb = 500 cm s~^ as the boundary condition at the 



outer (in radius) boundary rmax and assumed the 
radial velocity to be zero. A grid of 200 computed 
zones in radius was used in the numerical integration. 
An expansion in Legendre polynomials up to the 
number /max = 20 was used to solve the Poisson 
equation. 

A similar flow pattern was observed in all of these 
cases. Figures 1—3 show plots of the radial veloc- 
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ity against the radius (profiles) at consecutive times. 
As we see from ttiese figures, ttie matter inside ttie 
computed region in ail calculations falls through the 
transparent inner boundary to the center under the 
effect of gravity in a regime resembling free fall. The 
neutrino signal, which causes the mass of the PNS 
located at the center to decrease, manifests itself only 
in the existence of a modest maximum in the radial- 
velocity profile. In calculations [2, 3] with more mas- 
sive stars (2OM0 and 25Mq), the velocity maximum 
emerges immediately ahead of the front of a strong 
rarefaction wave: the maximum radial velocity is pos- 
itive but low; it does not exceed ~2.8 x 10'' cm s^^. 
In calculation [1] with a low-mass star (IIM0), the 
maximum lies in the range of negative velocities, so 
there are no positive radial velocities along the en- 
tire profile. We also see from Figs. 1—3 that in each 
case, a zone of positive velocities appears near the 
outer boundary of the computed region. This fact can 
be explained in terms of the boundary condition at 
the outer radius, which simulates a vacuum. This 
boundary condition has no fundamental effect on the 
properties of the flow as a whole; it only gives rise 
to a weak rarefaction wave. The main conclusion 
drawn from all calculations [1]— [3] is that almost all 
of the matter inside the computed region is accreted 
onto the central PNS. The plots of the radial velocity 
against the mass coordinate (profiles) show that the 
point of zero velocity (the point of separation) moves 
along the mass coordinate outward during the entire 
calculation. This implies that if the outer stellar en- 
velope is actually ejected, then the point of separation 
cannot be closer to the center than the inner boundary 
of the helium shell. 

The second series of our calculations was devoted 
to the effect of hot nucleon bubbles, which are most 
likely formed in the PNS corona, on the behavior 
of the onionlike structure of a massive star after the 
onset of gravitational iron-core collapse. In the initial 
data of these calculations, a hot bubble was specified 
in the following artificial way: the temperature in 
several adjacent cells near the inner boundary of 
the computed region was raised compared to the 
surrounding matter to the extent that, according to 
the adopted equation of state, the mass fractions 
of iron and helium nuclei became zero. The distri- 
bution of the remaining thermodynamic quantities 
was constructed to satisfy the hydrostatic equilibrium 
conditions ( 1 1 ), ( 12) as before. All calculations were 
performed with the same model of a 25Mq star; only 
the mass of the hot nucleon gas in the bubble was 
varied. Since the problem was solved in the one- 
dimensional approximation, the region of hot nucleon 
gas was actually a spherical shell rather than a bub- 
ble. The system of hydrodynamic equations (6)— (10) 
was integrated with the following initial data: [4] the 



minimum bubble radius = 8.901 x 10^ cm, the 
maximum bubble radius r^^x = 1-042 x 10^ cm, 
the hot-gas mass Mnb = 0.396 x IO'^Mq, and the 
mass of the matter inside the computed region Mjn = 
5.779M0; [5] r^? = 8.901 X 10^ cm, r^l^ = 1.248 x 
10^ cm, Mnb = 0.693 x W'^Mq, Mm = S.SSSMq; 
[6] r^B^ = 8.901 X lO'^ cm, r^^^ = 1.499 x 10^ cm, 
Mnb = 0.996 x lO'^M©, and Mm = 6.O36M0. Cal- 
culation [7] differs from case [6] only in that a positive 
velocity, Fnb = 1-734 x 10^ cm s~^, equal in mag- 
nitude to the local free-fall velocity was assigned to 
the matter in the cells that fell within the region of 
hot gas. In calculations [4]— [7], the parameters of 
the computed region were taken to be identical: the 
total stellar mass Mtot = 25M0, the iron-core mass 
Mcore = I-62OM0, the gravitational iron-core mass 
defect AMq = 4.386 x 10^^ g, the inner radius of the 
computed region rmin = 8.764 x lO''' cm, its outer 
radius rmax = 6.607 x 10^ cm, and the mass of the 
central region Mcenter = 1-OOO1M0. Thus, all of the 
latter parameters are identical to those in case [3]. 

In Figs. 4—7, the radial velocity is plotted against 
radius at consecutive times for cases [4]— [7], respec- 
tively. As we see from these figures, the shock wave 
produced through energy release in the nucleon bub- 
ble due to proton and neutron recombination into iron 
nuclei has an amplitude that is qualitatively propor- 
tional to the initial mass of the hot gas. In case [4], 
the shock wave is rather weak and the velocity of the 
post-shock matter barely reaches ~3 x 10^ cm s~^. 
In this case, a strong rarefaction wave is formed at the 
inner boundary of the computed region, as in the cal- 
culations that only took into account the gravitational 
effect of neutrino radiation, and the matter begins 
to be accreted onto the gravitating center at a high 
rate. We see from Fig. 8 that almost all of the matter 
that was initially inside the computed region passes 
through the inner boundary, whereas the escape of the 
matter through the outer boundary is negligible and is 
entirely attributable to the artificial rarefaction wave 
produced by a fixed boundary condition. The time it 
took for the shock wave to reach the helium shell is 
~13 s. In cases [6]— [7], the nucleon-gas mass is at 
a maximum (~O.O1M0). The appearance of a strong 
shock wave and envelope ejection were shown in both 
cases: the post-shock radial velocity is ~1.1 x 10^ 
and --^1.5 x 10^ cm s~^ for calculations [6] and [7], 
respectively. Therefore, the positive velocity (~1.7 x 
10^ cm s~^) imparted to the hot gas of the nucleon 
bubble at the initial time in case [7] has a marginal 
effect on the envelope ejection. The rarefaction wave 
formed at the inner computed boundary is found to 
be much weaker than that in calculation [4]. As a 
result, the amount of matter accreted onto the center 
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in these two calculations is modest: ~O.O45M0 and 
~O.O2OM0 for a zero and nonzero velocity of the 
hot gas, respectively. Almost all of the matter inside 
the computed region is ejected outward by the shock 
wave, as clearly confirmed by the plots in Figs. 9 



and 10. In case [6], it took ~6 s for the shock wave 
to reach the helium shell and ^^8 s for all of the matter 
to pass through the outer boundary; in case [7], both 

processes take about 4 s. 

The calculation with the initial data [5] is interme- 
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diate between case [4], on the one hand, and cases [6] 
and [7], on the other hand. A moderate rarefaction 
wave is formed near the inner boundary and the post- 
shoct; gas velocity is ~5 x 10*^ cm s~^ (Fig. 5). The 
shoclt wave reaches the helium shell in ~10 s. In this 



case, mass accretion at the inner computed boundary 
begins with an appreciable delay of ~12 s and the 
matter flows outward from the center at a rate that is 
approximately twice the rate of the gas accretion onto 
the PNS(Fig. 11). 
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Fig. 8. Total mass fluxes through the outer and inner 
boundaries of the computed region versus time for the 
calculations with data set [4]; (Mm is the total mass of 
the matter inside the computed region, Mdown is the total 
mass flux through the inner boundary, Mup is the total 
mass flux through the outer boundary, and Mtot = Mn + 
Miom + Mup is the total mass of the matter. 





Fig. 10. Same as Fig. 8 for data set [7]. 




Fig. 9. Same as Fig. 8 for data set [6]. 



Fig. 1 1 . Same as Fig. 8 for data set [5]. 
tend = 14.23 s; [7]: ^kin out = 5.46 x lO^i erg and 

tend = 11.03 S. 



In cases [4]— [7], nucleons are recombined into iron 
nuclei rapidly and in approximately the same time 
for all cases of initial data, which does not exceed 
~ 1 s. For cases [4]— [7], we also calculated the kinetic 
energy flux through the outer boundary in radius rmax- 
The following values were obtained for the total flux, 
which may prove to be of use in assessing the pos- 
sibility of emersion of the outer stellar layers: [4] the 
total kinetic energy flux through the outer boundary 
by the completion of the calculation E'kin out = 1-80 x 
10^^ erg and the completion time of the calculation 
tend = 16.84 s; [5]: ^kin out = 3.87 x lO^o erg and 
tend = 19.65 s; [6]: Ei^in^ui = 2.23 x lO^i erg and 



CONCLUSIONS 

Of course, the conclusions regarding the ejec- 
tion of the entire outer envelope of a massive star 
during the collapse of its iron core should have 
been reinforced by hydrodynamic calculations of the 
propagation of a strong shock wave through the 
helium shell and through the outer hydrogen— helium 
shell, certainly if they were preserved during the stellar 
evolution of this star. However, the velocity of the 
post-shock matter derived in the calculations with 
supercritical nucleon-bubble masses, which exceeds 
10^ cm s"-*^ for a radius of ~6.6 x 10^ cm, is, in 
turn, appreciably higher than the characteristic free- 
fall velocity with an inner mass of 6.8Mq (M = 
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Min + Mcenter — AMg), by more than a factor of 2 in 
case [6] (vff = f^^M.\ ' = 5.23 x 10^ cm/s J . 

The total outward kinetic energy flux through the 
outer boundary (see Sect. "Numerical Results") and 
the estimated recombination energy of the nucleon 
bubble simulated in the initial conditions of our 
problem suggest that complete ejection is possible. 
For the nucleon-bubble mass, we obtain Mnb = 

8 79 MeV 

0.996 X IQ-^Mq, so = Mnb x = 

mo 

1.68 X 10^° erg. This energy is close in absolute value 
to the gravitational energy of the outer stellar layers, 
~2.5 X 10^° erg. 

It should be particularly emphasized that our 
calculations yielded the critical nucleon-bubble mass 
(MNB)crit = ^NB ([5]) = 0.7 X W'^ Mq (scc casc 
[5]). It should be remembered that here, we deal with 
the most massive star with Mms = 25Mq. For less 
massive stars, (Af^B)crit most likely much lower, al- 
though this suggestion needs to be verified by further 
calculations. Can such nucleon bubbles, which assist 
in the envelope ejection more than the other gravita- 
tional effect of energy losses in the form of a neutrino 
signal taken into account in our calculations, be 
formed? In the one-dimensional hydrodynamic theory 
of gravitational collapse (in many studies with an 
allowance made for neutrino radiation, including that 
at the neutrinosphere formation stage, which virtually 
coincides with the second collapse stage mentioned 
in the Introduction), a rather thick spherical shell has 
long been shown to appear above the neutrinosphere, 
which is semitransparent to neutrino radiation [see 
Imshennik and Nadyozhin (1982) for a review]. Bur- 
rows and Goshy (1993) showed that this shell has 
a quasi-steady-state structure in the typical case of 
a frozen accretion wave whose front represents the 
outer shell boundary. It is important that this shell is 
composed of a Boltzmann gas of free nucleons and 
electron— positron plasma. Imshennik (2002) refined 
such quasi-steady-state solutions and showed its 
strong convective instability virtually in the entire 
solution region that encompasses all the possible 
conditions for gravitational collapse or, to be more 
precise, its second stage (post-shock accretion). The 
typical masses in these shells called neutrino PNS 
coronas are 1O~''^M0 < Mcpns < 1O~^M0; i.e., they 
are roughly equal to the critical mass M^b derived 
above. The possibility of such nucleon bubbles being 
ejected due to the nonlinear development of convec- 
tion is supported by the well-known considerations 
given by Burrows et al. (1994) based on their calcu- 
lations of a two-dimensional model for gravitational 
collapse and by the interesting self-similar solutions 



of a neutrino wind that drags nucleon bunches into 
the outer stellar envelopes (Thompson et al. 2001). 
Here, we demonstrated the hydrodynamic effect of the 
ejection of the outer stellar envelope in the rough one- 
dimensional approximation. However, to reach the 
final conclusion regarding the ejection of the entire 
outer envelope, strictly speaking, requires further cal- 
culations of the passage of the shock wave obtained 
in our calculations through the helium shell and the 
outer hydrogen— helium shell. Still, it is qualitatively 
clear that a significant fraction of the kinetic energy 
accumulated behind the shock front will be spent on 
overcoming the gravitational binding energy of these 
outermost layers. For this reason, the resulting mean 
ejection velocities of the entire envelope can be much 
lower than the calculated velocities; i.e., the envelope 
will be essentially emersed rather than ejected. 
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